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Abstract 

Background: An efficient and reliable parameter estimation method is essential for the creation of biological 
models using ordinary differential equation (ODE). Most of the existing estimation methods involve finding the 
global minimum of data fitting residuals over the entire parameter space simultaneously. Unfortunately, the 
associated computational requirement often becomes prohibitively high due to the large number of parameters 
and the lack of complete parameter identifiability (i.e. not all parameters can be uniquely identified). 

Results: In this work, an incremental approach was applied to the parameter estimation of ODE models from 
concentration time profiles. Particularly, the method was developed to address a commonly encountered 
circumstance in the modeling of metabolic networks, where the number of metabolic fluxes (reaction rates) 
exceeds that of metabolites (chemical species). Here, the minimization of model residuals was performed over a 
subset of the parameter space that is associated with the degrees of freedom in the dynamic flux estimation from 
the concentration time-slopes. The efficacy of this method was demonstrated using two generalized mass action 
(GMA) models, where the method significantly outperformed single-step estimations. In addition, an extension of 
the estimation method to handle missing data is also presented. 

Conclusions: The proposed incremental estimation method is able to tackle the issue on the lack of complete 
parameter identifiability and to significantly reduce the computational efforts in estimating model parameters, 
which will facilitate kinetic modeling of genome-scale cellular metabolism in the future. 
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Background 

The estimation of unknown kinetic parameters from 
time-series measurements of biological molecules is a 
major bottleneck in the ODE model building process in 
systems biology and metabolic engineering [1]. The ma- 
jority of current estimation methods involve simultan- 
eous (single-step) parameter identification, where model 
prediction errors are minimized over the entire param- 
eter space. These methods often rely on global optimization 
methods, such as simulated annealing, genetic algorithms 
and other evolutionary approaches [1-3]. The problem of 
obtaining the best-fit parameter estimates however, is typic- 
ally ill-posed due to issues related with data informative- 
ness, problem formulation and parameter correlation, all of 
which contribute to the lack of complete parameter 
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identifiability. Not to mention, finding the global minimum 
of model residuals over highly multidimensional parameter 
space is challenging and can become prohibitively expen- 
sive to perform on a computer workstation, even for tens 
of parameters. 

Here, we consider the modeling of cellular metabolism 
using the canonical power-law formalism, specifically the 
generalized mass action (GMA) systems [4,5]. The power- 
law formalism has many advantages, which have been 
detailed elsewhere [1,6]. Notably, power laws have a rela- 
tively simple structure that permits algebraic manipulation 
in the logarithmic scale, but nonetheless is capable of de- 
scribing essentially any nonlinearity. Regulatory interac- 
tions among metabolites can also be described 
straightforwardly through the kinetic order parameters, 
establishing an equivalence between structural identifica- 
tion and parametric estimation. However, the number of 
parameters increases proportionally with the number of 
metabolites and fluxes, leading to a large-scale parameter 
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identification problem, one where single-step estimation 
methods often struggle to converge. 

The integration of ODE often constitutes a major part 
of the computational cost in the parameter estimation, 
especially when the ODE model is stiff [7]. While stiff- 
ness can genuinely arise due to a large time scale separ- 
ation of the reaction kinetics in the real system, stiff 
ODEs could also result from unrealistic combinations of 
parameter values during the parameter optimization 
procedure, especially when a global optimizer is used. 
The parameter estimation of ODE models using power- 
law kinetics is particularly prone to stiffness problem 
since many of the unknown parameters are the expo- 
nents of the concentrations. For this reason, alternative 
formulations have been proposed that avoid these ODE 
integrations either completely [7,8] or partially [9-11]. 
Particularly, computational cost could be significantly 
reduced by decomposing the estimation problem into 
two phases, starting with the calculation of dynamic re- 
action rates or fluxes from the slopes of concentration 
data, followed by the least square regressions of kinetic 
parameters [12-14]. In this case, the final parameter esti- 
mation is done one flux at a time, each involving only a 
handful of parameters and thus, the global minimum so- 
lution can be either computed analytically (for example, 
when using log-linear power-law flux functions) or 
determined efficiently. Moreover, as the first estimation 
phase (flux estimation) depends only on the assumption 
of the topology of the metabolic network, the flux esti- 
mates can subsequently be used to guide the selection of 
the most appropriate flux functions for the second phase 
or to detect inconsistencies in the assumed topology of 
the network separately from the flux equations [14]. 
However, the application of this method requires the 
number of metabolites to be equal to or larger than that 
of fluxes, so that the flux estimation can result in a 
unique solution. Since the reverse situation is more 
commonly encountered in the typical metabolic net- 
works, a generalization of this incremental estimation 
approach becomes the main focus in this study. 

As noted above, the new parameter estimation method 
in this work is built on the concept of incremental iden- 
tification [12,13] or dynamical flux estimation (DFE) 
method [14,15]. The proposed method provides two new 
contributions: (1) an ability to handle the more general 
scenario, where the number of reactions exceeds that of 
the metabolites and (2) high numerical efficiency 
through the reduction of the parameter search space. 
Specifically, two parameter estimation formulations are 
proposed with objective functions that depend on model 
prediction errors of metabolite concentrations and of 
concentration time-slopes. An extension of this strategy 
to circumstances where concentration data of some 
metabolites are missing is also presented. The proposed 



method is applied to two previously published GMA 
models and compared with single-step estimation meth- 
ods, in order to demonstrate its efficacy. 

Methods 

The generalized mass action model of cellular metabol- 
ism describes the mass balance of metabolites, taking 
into account all metabolic influxes and effluxes and their 
stoichiometric ratios, as follows: 

dX(t, p)/dt = X(f, p) = Sv(X, p), (1) 

where X(£,p) is the vector of metabolic concentration 
time profiles, S £ R m x n is the stoichiometric matrix for 
m metabolites that participate in n reactions, and v(X,p) 
denotes the vector of metabolic fluxes (i.e. reaction 
rates). Here, each flux is described by a power-law equa- 
tion: 

V/(X,p) = r /Qxf\ (2) 

i 

where jj is the rate constant of the y-th flux and is the 
kinetic order parameter, representing the influence of 
metabolite X t on the y-th flux (positive: X t is an activating 
factor or a substrate, negative: X t is an inhibiting factor). In 
incremental parameter identification, a data pre-processing 
step (e.g. smoothing or filtering) is usually applied to the 
noisy time-course concentration data X m fe), in order to 
improve the time-slope estimates X m (tk) . Subsequently, 
the dynamic metabolic fluxes vfe) are estimated from 
Equation (1) by substituting X(t) with X m (tj < ). Finally, the 
kinetic parameters associated with the y-th flux (i.e. jj an d 
fji's) can be calculated using a least square regression of 
the power law flux function in Equation (2) against 
the estimated Vj(t k ). Note that for GMA models, the 
least square parameter regressions in the last step are 
linear in the logarithmic scale and thus, can be 
performed very efficiently. 

A unique set of dynamic flux values vfe) can only be 
computed from X m (tk) = Sv(^), when the number of 
metabolites exceeds that of fluxes. However, a metabol- 
ite in general can participate in more than one meta- 
bolic flux (m<n). In such a situation, there exist an 
infinite number of dynamic flux combinations v(t k ) that 
satisfy X m (tk) = Sv(fe). The dimensionality of the set of 
flux solutions is equal to the degree of freedom (DOF), 
given by the difference between the number of fluxes 
and the number of metabolites: n DO F= n-m >0 (assum- 
ing S has a full row rank, i.e. there is no redundant 
ODE in Equation (1)). The positive DOF means that 
the values of n DOF selected fluxes can be independently 
set, from which the remaining fluxes can be computed. 
This relationship forms the basis of the proposed 
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estimation method, in which the model goodness of fit 
to data is optimized by adjusting only a subset of para- 
meters associated with the independent fluxes above. 

Specifically, we start by decomposing the fluxes into 
two groups: vfe) = [ v 7 (^) r \D(fk) T ] T > where the 
subscripts / and D denote the independent and 
dependent subset, respectively. Then, the parameter vec- 
tor p and the stoichiometric matrix S can be structured 
correspondingly as p = [ p 7 p^ ] and S = [ S 7 S D ]. The rela- 
tionship between the independent and dependent fluxes 
can be formulated by rearranging X w (^) = Sv(fe) into: 



vl>(£/c) = S D 1 [x w (£/ c ) - S/V/(X w (fe),p 7 )j . 



(3) 



In this case, given p 7 , one can compute the independent 
fluxes v 7 (X m (£/ c ),p 7 ) using the concentration data X m (t/X 
and subsequently obtain Voitk) from Equation (3). 
Finally, p^ can be estimated by a simple least square 
fitting of V£>(X m (£/ c ),p£>) to the computed V£>(£/ c ) one flux 
at a time, when there are more time points than the 
number of parameters in each flux. 

In this study, two formulations of the parameter esti- 
mation of ODE models in Equation (1) are investigated, 
involving the minimization of concentration and slope 
errors. The objective function for the concentration 
error is given by 



<Mp,X) = 



\ 



i K 

—Y, PUfc) - x(t k , P )] T [x m (t k ) - x(t k , P ) 



(4) 



and that for the slope error is given by 



Os(p,X) 



1 r V 

k=l 

\[± m (t k ) -Sv(X OT (^),p)], 



(5) 



where K denotes the total number of measurement time 
points and X(^,p) is the concentration prediction (i.e. 
the solution to the ODE model in Equation (1)). Figure 1 
describes the formulation of the incremental parameter 
estimation and the procedure for computing the object- 
ive functions. Note that the computation of O c requires 
an integration of the ODE model and thus, the estima- 
tion using this objective function is expected to be com- 
putationally costlier than that using O s . On the other 
hand, metabolic mass balance is only approximately 
satisfied at discrete time points during the parameter 
estimation using O s , as the ODE model is not integrated. 



minO(p,x) 

\p D e[L,U] 



s.t. 



1. Calculate independent fluxes 



rjTlx? ji 



(n-m)xl 



2. Solve for dependent fluxes j 

x„ft K s, sj££>] | 

V,) = S 0 - , [*.((,)-S,t,((,)]| 




4. Compute objective function 



3. Estimate the parameters associated 
with dependent fluxes 



lnv D (t k ) = 



ln r .+I/„.lnX i (0 



Figure 1 Flowchart of the incremental parameter estimation. 
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There are several important practical considerations in 
the implementation of the proposed method. The first 
consideration is on the selection of the independent 
fluxes. Here, the set of these fluxes is selected such that 
(i) the mxm submatrix is invertible, (ii) the total 
number of the independent parameters p 7 is small, and 
(iii) the prior knowledge of the corresponding p 7 is maxi- 
mized. The last two aspects should lead to a reduction 
in the parameter search space and the cost of finding 
the global optimal solution of the minimization problem 
in Figure 1. The second consideration is regarding con- 
straints in the parameter estimation. Biologically relevant 
values of parameters are often available, providing lower 
and/or upper bounds for the parameter estimates. In 
addition, enzymatic reactions in the ODE model are often 
assumed to be irreversible and thus, dynamic flux esti- 
mates are constrained to be positive. Hence, the parameter 
estimation involves a constrained minimization problem, 
for which many global optimization algorithms exist. 

So far, we have assumed that the time-course concen- 
tration data are available for all metabolites. However, 
the method above can be modified to accommodate 
more general circumstances, in which data for one or 
several metabolites are missing. In this case, the ODE 



model is first rewritten to separate the mass balances 
associated with measured and unmeasured metabolites, 
such that 



(6) 



where the subscripts M and U refer to components that 
correspond to measured and unmeasured metabolites, 
respectively. Again, if the fluxes are split into two cat- 
egories v 7 and V£, as above, the following relationship still 
applies for the measured metabolites: 



X(f,p) = 




(t,p) = 


Sat 


v(X M ,X u ,p) 








.V 





l(tk)] 



(7) 



Naturally, the degree of freedom associated with the 
dynamic flux estimation is higher by the number of 
component in than before. Figure 2 presents a 
modification of the parameter estimation procedure in 
Figure 1 to handle the case of missing data, in which an 
additional step involving the simulation of unmeasured 
metabolites X^ = Suv(X M , X^, p) will be performed. In 
this integration, X M is set as an external variable, whose 
time-profiles are interpolated from the measured 



minO(p,x) 

X v >0 



s.tJ 



v D (^)>0 
P D e[L,C/] 



j 1. Simulate unmeasured metabolites j 
\ 2. Calculate independent fluxes j j 5. Compute objective function 



(n-m)xl 



I 



4. Estimate the parameters associated 
with dependent fluxes 



3. Solve for dependent fluxes 

| X m (4) = [S;,M S DM ] 



lnv D (f,) = 



v,(4) 

v D (fJ = S D>M -'[X m (f k )-S ; , M v ; (fj] 



ln r .+I/. ; lnX(g 



Figure 2 Flowchart of the incremental parameter estimation when metabolites are not completely measured. 
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concentrations. The set of independent fluxes v 7 are now 
selected to include all fluxes that appear in and those 
that lead to a full column ranked S D>M . If S D>M is a non- 
square matrix, then a pseudo-inverse will be done in 
Equation (7). Of course, the same considerations men- 
tioned above are equally relevant in this case. Note that 
the initial conditions of X u will also need to be estimated. 



Results 

Two case studies: a generic branched pathway [7] and 
the glycolytic pathway of L. lactis [16], were used to 
evaluate the performance of the proposed estimation 
method. In addition, simultaneous estimation methods 
employing the same objective functions in Equations (4) 
and (5) were applied to these case studies, to gauge the 
reduction in the computational cost from using the pro- 
posed strategy. In order to alleviate the ODE stiffness 
issue, parameter combinations that lead to a violation in 
the MATLAB (ode 15s) integration time step criterion is 
assigned a large error value (O c = 10 3 for the branched 
pathway and 10 5 for the glycolytic pathway). Alterna- 
tively, one could also set a maximum allowable integra- 
tion time and penalize the associated parameter values 
upon violation, as described above. In this study, the 
optimization problems were solved in MATLAB using 
publicly available eSSM GO (Enhanced Scatter Search 
Method for Global Optimization) toolbox, a population- 



based metaheuristic global optimization method incorp- 
orating probabilistic and deterministic strategies [17,18]. 
The MATLAB codes of the case studies below are avail- 
able in Additional file 1. Each parameter estimation was 
repeated five times to ensure the reliability of the global 
optimal solution. Unless noted differently, the iterations 
in the optimization algorithm were terminated when the 
values of objective functions improve by less than 0.01% 
or the runtime has exceeded the maximum duration (5 
days). 

A generic branched pathway 

The generic branched pathway in this example consists 
of four metabolites and six fluxes, describing the trans- 
formations among the metabolites (double-line arrows), 
with feedback activation and inhibition (dashed arrows 
with plus or minus signs, respectively), as shown in 
Figure 3 A. The GMA model of this pathway is given in 
Figure 3B, containing a total of thirteen rate constants 
and kinetic orders. This model with the parameter 
values and initial conditions reported previously [7] were 
used to generate noise-free and noisy time-course con- 
centration data (i.i.d additive noise from a Gaussian dis- 
tribution with 10% coefficient of variation). The noisy 
data were smoothened using a 6-th order polynomial, 
which provided the best relative goodness of fit among 
polynomials according to Akaike Information Criterion 
(AIC) [19] and adjusted R 2 [20]. Subsequently, time- 



B 








"i 


K 




0 






0 






0 



0 0 



0 

-1 
1 

0 









v l 


0 


-1 


0" 


V 2 


0 


0 


0 


V 3 


-1 


0 


0 


V 4 


0 


1 


-1 


V 5 








_ V 6 



v, = y x x 0 x; f » 
v 2 =r 2 x/» 

v 3 = r3 X/» 
v 4 = y A X z u X/« 
v 5 = y 5 X/» 

n = r 6 x/« 

v X 0 = 0.6 

Figure 3 A generic branched pathway. (A) Metabolic pathway map and (B) the GMA model equations [7]. 



X x (0 = 1-40 
X 2 (f 0 ) = 2.70 
X 3 (0 = 1-20 
X 4 (0 = 0.40 



Jia et al. BMC Systems Biology 201 2, 6:142 
http://www.biomedcentral.eom/1 752-0509/6/1 42 



Page 6 of 12 



slopes of noise-free and smoothened noisy data were com- 
puted using the central finite difference approximation. 

Here, v 2 and v 6 were chosen as the independent fluxes 
as they comprise the least number of kinetic parameters 
and lead to an invertible S^,. The two rate constants and 
two kinetic orders were constrained to within [0,25] and 
[0,2], respectively. In addition, all the reactions are 
assumed to be irreversible. 

Table 1 compares simultaneous and incremental par- 
ameter estimation runs using noise-free data, employing 
the two objective functions above. Regardless of the ob- 
jective function, the proposed incremental approach sig- 
nificantly outperformed the simultaneous estimation. 
When using the concentration-error minimization, sim- 
ultaneous optimization met great difficulty to converge 
due to stiff ODE integrations. Only one out of five 
repeated runs could complete after relaxing the conver- 
gence criteria of the objective function to 1%, while the 
others were prematurely terminated after the prescribed 
maximum runtime of 5 days. In contrast, the proposed 
incremental estimation was able to find a minima of O c 
in less than 96 seconds on average with good concentra- 
tion fit and parameter accuracy (see Figure 4A and 
Table 1). By avoiding ODE integrations using O s , the 
simultaneous estimation of parameters could be com- 
pleted in roughly 10 minutes duration, but this was 
much slower than the incremental estimation using O c . 
In this case, the incremental method was able to con- 
verge in below 2 seconds or over 250 times faster. The 
goodness of fit to concentration data and the accuracy 
of parameter estimates were relatively equal for all three 
completed estimations (see Figure 4B and Table 1). The 
parameter inaccuracy in this case was mainly due to the 
polynomial smoothing of the concentration data, since 
the same estimations using the analytical values of the 
slopes (by evaluating the right hand side of the ODE 
model in Equation (1)) could give accurate parameter 
estimates (see Additional file 2: Table SI). 
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Figure 4 Simultaneous and incremental estimation of the 
branched pathway using in silico noise-free data (x). (A) 

concentration predictions using parameter estimates from 
incremental method by 0 C minimization ( — ); (B) concentration 
predictions using parameter estimates from simultaneous method 
(o) and proposed method ( — ) by 0 S minimization. 



Table 2 provides the results of the same estimation 
procedures as above using noisy data. Data noise led to 
a loss of information and an expected decline in the par- 
ameter accuracy. Like before, the simultaneous estima- 
tion using ® c met stiffness problem and three out of 
five runs did not finish within the five-day time limit. 
The incremental approach using either one of the ob- 
jective functions offered a significant reduction in the 
computational time over the simultaneous estimation 
using O s , while providing comparable parameter accur- 
acy and concentration and slope fit (see Figure 5 and 



Table 1 Parameter estimations of the branched pathway model using noise-free data 






Simultaneous method 


Incremental method 




min 0£ 


min (Dg 


min O c 


min O s 


CPU time (sec) a 


56.00 h 


620.81 ±64.30 


95.95 ±11.09 


1.56 ±0.19 


eSSM GO iterations 


323 


4390 + 391 


14±4 


10 + 2 


Parameter error (%) 


49.10 


36.91% ±1.09 


21.56%±7.57x10" 2 


36.85% ± 6.48 x 10" 3 


< 


4.54 x 10~ 3 


6.54 x 10" 3 ±5.20x 10~ 5 


4.03 X10" 3 ±6.22x10" 8 


6.00 x 10" 3 ±5.05x 10~ 7 


CDs 


7.01 x 10" 2 


2.72 x 10" 2 ± 1.09 x 10" 5 


3.92 x 10" 2 ± 9.86 X10" 6 


2.76 x 10" 2 ±4.46x 10~ 10 



a. CPU time was based on a workstation with dual Intel Quad-Core 2.83 GHz processors. 

b. Only one out of five runs completed with a relative improvement of the objective function below 1% between iterations. The rest did not converge within the 
5-day time limit after iterating for 583, 989, 777, and 661 times. The corresponding CD C at termination were 4.85x 10~ 2 , 1.39 x 10~ 2 , 1.75 x 10~ 2 and 3.75 x 10~ 2 , 
respectively. 

c. Mean value ± standard deviation out of five repeats. 

d. Root mean square error of model predictions, where the underlined value refers to the objective function of the minimization. 
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Table 2 Parameter estimations of the branched pathway model using noisy data 

Simultaneous method 



Incremental method 



min Of 



min O s 



min O c 



min O s 



CPU time (sec) 


1 7.86 h 


534.83 ±22.1 2 


71.88 ±6.33 


1.17 + 0.12 




44.63 h 








eSSM GO iterations 


254 


3494 ±348 


12 + 2 


10 + 3 




426 








Parameter error (%) 


75.42 


54.36 ± 4.47 


75.77 ±6.11 x 10" 3 


51.15 + 1.38X 10" 3 




34.98 








0c 


3.62 X10" 2 


6.06x10" 2 ±1.14x10" 3 


3.52 x 10" 2 ± 9.50 X10" 9 


4.76 x 10" 2 ± 3.81 x 10~ 7 




3.27 X10" 2 








0s 


2.06 x 10" 1 


1.34x10" 1 ±6.02x10" 4 


1.64X 10" 1 ±2.23x 10" 5 


1.38X 10~ 1 ±2.26x 10" 10 




1.60X10" 1 









a. Two out of five runs completed with a relative improvement of the objective function below 1% between iterations. The rest did not converge within the 5-day 
time limit after iterating for 805, 699, and 568 times. The corresponding CD C at termination were 4.08 x 10~ 2 , 5.05 x 10~ 2 and 6.25 x 10~ 2 , respectively. 



Table 2). In this example, data noise did not affect the 
computational cost in obtaining the (global) minimum 
of the objective functions. 

Finally, the estimation strategy described in Figure 2 
was applied to this example using noise-free data and as- 
suming X 3 data were missing. Fluxes v 3 and v 4 that ap- 
pear in X 3 were chosen to be among the independent 
fluxes and flux Vj was also added to the set such that the 
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Figure 5 Simultaneous and incremental estimation of the 
branched pathway using in silico noisy data (x). (A) 

concentration predictions using parameter estimates from 
incremental method by 0 C minimization ( — ); (B) concentration 
predictions using parameter estimates from simultaneous method 
(o) and proposed method ( — ) by O s minimization. 



dependent fluxes can be uniquely determined from 
Equation (7). In addition to the parameters associated 
with the aforementioned fluxes, the initial condition X 3 (t 0 ) 
was also estimated. The bounds for the rate constants and 
kinetic orders were kept the same as above, while the 
initial concentration was bounded within [0, 5]. 

Table 3 summarizes the parameter estimation results. 
Four out of five repeated runs of O c simultaneous 
optimization were again prematurely terminated after 5 
days. Meanwhile, the rest of the estimations could pro- 
vide reasonably good data fitting with the exception of 
fitting to X 3 data as expected (see Figure 6). Like data 
noise, missing data led to increased inaccuracy of the 
parameter estimates, regardless of the estimation meth- 
ods. Finally, the computational speedup by using the 
incremental over the simultaneous estimation was sig- 
nificant, but was lower than in the previous runs due to 
the additional integration of and the larger number 
of independent parameters. The detailed values of the 
parameter estimates in this case study can be found in 
the Additional file 2: Tables S2 and S3. 

The glycolytic pathway in Lactococcus. lactis 

The second case study was taken from the GMA model- 
ing of the glycolytic pathway in L. lactis [16], involving six 
internal metabolites: glucose 6-phosphate (G6P) - X h 
fructose 1, 6-biphosphate (FBP) - X 2 , 3-phosphoglycerate 
(3-PGA) - X 3) phosphoenolpyruvate (PEP) - X 4 , Pyruvate 
- X 5 , Lactate - X 6 , and nine metabolic fluxes. In addition, 
external glucose (Glu), ATP and Pi are treated as off-line 
variables, whose values were interpolated from measure- 
ment data. The pathway connectivity is given in Figure 7A, 
while the model equations are provided in Figure 7B. 

The time-course concentration dataset of all metabo- 
lites were measured using in vivo NMR [21,22], and 
smoothened data used for the parameter estimations 
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Table 3 Parameter estimations of the branched pathway model using noise-free data with X 3 missing 

Simultaneous method Incremental method 



min Oc min O s min O c min O s 



CPU time (sec) 


85.03 h 


4002.01 ±69 


6.11 


1 404.22 ± 120.71 


445.47 ± 35.94 


eSSM GO iterations 


308 


365 ±91 




67 ± 10 


48 ±10 


Parameter error (%) 


71.90 


43.50 ± 2.34 




68.85 ±4.57 


40.47 ± 0.59 


0c 


4.54 x 10" 3 


6.46x10" 3 ± 


4.08 X10" 4 


3.38x10" 3 ±1.14x10" 4 


5.94 x 10" 3 ± 3.23 X10" 5 


0s 


1.03 


2.99x10" 2 ± 


3.82 X10" 4 


8.32 x 10" 2 ± 4.04 X10" 3 


2.94 x 10" 2 ± 2.77 X10" 6 



a. Only one out of five runs completed with a relative improvement of the objective function below 1% between iterations. The rest did not converge within the 
5-day time limit after iterating for 471, 435, 863 and 786 times. The corresponding (D c at termination were 4.99x 10" 2 , 4.92 x 10" 2 , 1.17 x 10~ 2 and 1.57 x 10" 2 , 
respectively. 



below were shown in Figure 8. The raw data has been 
filtered previously [16], and these smoothened data for 
all metabolites but X 6 , were directly used for the concen- 
tration slope calculation in this case study. In the case of 
X 6 , a saturating Hill- type equation: k 2 t n I (k 2 + t n ) where 
t is time and the constants kj, k 2 , n are smoothing para- 
meters, was fitted to the filtered data to remove unrealis- 
tic fluctuations. The central difference approximation 
was also adopted to obtain the time-slope data. 

Fluxes v 4) v 7 and v 9 were selected as the DOF, again to 
give the least number of p 7 and to ensure that is in- 
vertible. All rate constants were constrained to within 
[0, 50], while the independent and dependent kinetic 
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Figure 6 Simultaneous and incremental estimation of the 
branched pathway with missing X 3 : in silico noisy-free data (x). 

(A) concentration predictions using parameter estimates from 
incremental method by O c minimization ( — ); (B) concentration 
predictions using parameter estimates from simultaneous method 
(o) and proposed method ( — ) by O s minimization. 



orders were allowed within [0, 5] and [-5, 5], respectively. 
The difference between the bounds for the independent 
and dependent kinetic orders was done on purpose to 
simulate a scenario where the signs of the independent 
kinetic orders were known a priori. 

Table 4 reports the outcome of the single-step and in- 
cremental parameter estimation runs using O c and O s . 
The values of the parameter estimates are given in the 
Additional file 2: Table S4. Like in the previous case 
study, there was a significant reduction in the estimation 
runtime by using the proposed method over the simul- 
taneous estimation, with comparable goodness of fit in 
concentration and slope. None of the five repeats of O c 
simultaneous minimization converged within the five- 
day time limit, even after relaxing the convergence cri- 
teria of the objective function to 1%. On the other hand, 
the incremental estimation using O c was not only able 
to converge, but was also faster than the simultaneous 
estimation of O s that did not require any ODE integra- 
tion. The incremental estimation using O c was able to 
provide parameters with the best overall concentration 
fit (see Figure 8), despite having a large slope error. 
Finally, minimizing O s does not guarantee that the 
resulting ODE is numerically solvable, as was the case of 
simultaneous estimation, due to numerical stiffness. But 
the incremental parameter estimation from minimizing 
® s can produce solvable ODEs with good concentration 
and slope fits. 

Discussion 

In this study, an incremental strategy is used to develop 
a computationally efficient method for the parameter 
estimation of ODE models. Unlike most commonly used 
methods, where the parameter estimation is performed 
to minimize model residuals over the entire parameter 
space simultaneously, here the estimation is done in two 
incremental steps, involving the estimation of dynamic 
reaction rates or fluxes and flux-based parameter 
regressions. Importantly, the proposed strategy is 
designed to handle systems in which there exist extra 
degrees of freedom in the dynamic flux estimation, 
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_v 8 ^ (X 5 ) 

2 NAD+ 2 NADH 





= 5.73 




= 24.54 




= 7.50 


xM 


= 2.60 


X 5 {'o) 


= 12.72 


xA'o) 


= 0.43 



when the number of metabolic fluxes exceeds that of 
metabolites. The positive DOF means that there exist 
infinitely many solutions to the dynamic flux estimation, 
which is one of the factors underlying the parameter 
identifiability issues plaguing many estimation problems 
in systems biology [23,24]. 

The main premise of the new method is in recognizing 
that while many equivalent solutions exist for the dynamic 
flux estimation, the subsequent flux-based regression will 
give parameter values with different goodness -of-f it, as 
measured by O c or O s . In other words, given any 



two dynamic flux vectors v(t k ) satisfying X w (^) = 
Sv(^), the associated parameter pairs (p 7 , p^,) may not 
predict the slope or concentration data equally well, due 
to differences in the quality of parameter regression for 
each \(t k ). Also, because of the DOF, the minimization of 
model residuals needs to be done only over a subset of 
parameters that are associated with the flux degrees of 
freedom, resulting in much reduced parameter search 
space and correspondingly much faster convergence to 
the (global) optimal solution. The superior performance of 
the proposed method over simultaneous estimation was 
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Figure 8 Incremental estimation of the L lactis model: Experimental data (x) compared with model predictions using parameters from 
concentration error minimization ( — ) and slope error minimization ( — ). 
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Table 4 Parameter estimations of the L lactis model 


Simultaneous method 




Incremental method 




min Oc 


min O s 


min O c 


min O s 


CPU time (sec) >5 days 


3476.89 ± 349.63 


976.72 ±31.01 


20.82 ±2.71 


eSSM GO iterations — 


1662 ±282 


4±1 


33 ±7 


Oc 


Stiff ODE 


220 ±8.81 X10" 3 


6.1 8 ± 7.28 x 10~ 2 


Os 


2.67 ± 1.93 x 10" 4 


1.51 x10 3 ± 52.50 


5.79 ± 9.62 x 10" 4 



a. None of five runs finished with a relative improvement of the objective function below 1% within the 5-day time limit, after iterating for 60, 147, 93, 79 and 31 
times. The corresponding (D c at termination were 9.31, 7.57, 8.77, 9.39 and 12.9, respectively. 



convincingly demonstrated in the two GMA modeling 
case studies in the previous section. The minimization of 
slope error, also known as slope-estimation-decoupling 
strategy method [7], is arguably one of the most com- 
putationally efficient simultaneous methods. In this 
strategy, the parameter fitting essentially constitutes a 
zero-finding problem and the estimation can be done 
without having to integrate the ODEs. Yet, the incre- 
mental estimation could offer more than two orders 
of magnitude reduction in the computational time 
over this strategy. 

There are many factors, including data-related, model- 
related, computational and mathematical issues, which 
contribute to the difficulty in estimating kinetic para- 
meters of ODE models from time-course concentration 
data [1]. Each of these factors has been addressed to a 
certain degree by using the incremental identification 
strategy presented in this work. For example, in data- 
related issues, the proposed method can be modified to 
handle the absence of concentration data of some meta- 
bolites, as shown in Figure 2. Nevertheless, the method 
is neither able nor expected to resolve the lack of 
complete parameter identifiability due to insufficient 
(dynamical) information contained in the data [23,24]. 
As illustrated in the first case study, single-step and in- 
cremental approaches provided parameter estimates 
with similar accuracies, which expectedly deteriorated 
with noise contamination and loss of data. 

The appropriateness of using a particular mathemat- 
ical formulation, like power law, is an example of model- 
related issues. As discussed above, this issue can be 
addressed after the dynamic fluxes are estimated, where 
the chosen functional dependence of the fluxes on a spe- 
cific set of metabolite concentrations can be tested prior 
to the parameter regression [14]. Next, the com- 
putational issues associated with performing a global 
optimization over a large number of variables and the 
need to integrate ODEs have been mitigated in the pro- 
posed method by performing optimization only over the 
independent parameter subset and using a minimization 
of slope error, respectively. Finally, in this work, we have 
also addressed a mathematical issue related to the 
degrees of freedom that exist during the inference of 



dynamic fluxes from slopes of concentration data. How- 
ever, extra degrees of freedom (mathematical redundan- 
cies) are also expected to influence the second step of 
the method, i.e. one-flux-at-a-time parameter estimation. 
For (log)linear regression of parameters in GMA models, 
such redundancy will lead to a lack of full column rank 
of the matrix containing the logarithms of concentration 
data X m fe) and thus, can be straightforwardly detected. 

The proposed estimation method has several weak- 
nesses that are common among incremental estimation 
methods. As demonstrated in the first case study, the ac- 
curacy of the identified parameter relies on the ability to 
obtain good estimates of the concentration slopes. Direct 
slope estimation from the raw data, for example using 
central finite difference approximation, is usually not ad- 
visable due to high degree of noise in the typical bio- 
logical data. Hence, pre-smoothing of the time-course 
data is often required, as done in this study. Many algo- 
rithms are available for such purpose, from simplistic 
polynomial regression and splines to more advanced 
artificial neural network [7,25] and Whittaker-Eilers 
smoother [26,27]. If reliable concentration slope esti- 
mates are not available, but bounds for the slope values 
can be obtained, then one can use interval arithmetic to 
derive upper and lower limits for the dependent fluxes 
and parameters using Equation (3) (or Equation (7) [28]. 
When the objective function involves integrating the 
model, validated solution to ODE with interval para- 
meters can be used to produce the corresponding upper 
and lower bounds of concentration predictions [29]. 
Finally, the estimation can be reformulated, for example 
by minimizing the upper bound of the objective. 

In addition to the drawback discussed above, the pro- 
posed strategy requires a priori knowledge about the 
topology of the network. For cellular metabolism, such 
information has become more readily available as 
genome-scale metabolic network of many important 
organisms, including human, E. coli and S. cereviseae, 
have been and are continuously being reconstructed 
[30]. For other networks, many algorithms also exist for 
the estimation of network topology based on time-series 
concentration data, including Bayesian network infer- 
ence, transfer entropy, and Granger causality [31-33]. 
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Conclusions 

The estimation of kinetic parameters of ODE models 
from time-course concentration data remains a key 
bottleneck in model building in systems biology. The 
lack of complete parameter identifiability has been 
blamed as the root cause of the difficulty in such estima- 
tion. In this study, a new incremental estimation method 
is proposed that is able to overcome the existence of 
extra degrees of freedom in the dynamic flux estimation 
from concentration slopes and to significantly reduce 
the computational requirements in finding parameter 
estimates. The method can also be applied, after minor 
modifications, to circumstances where concentration 
data for a few molecules are missing. While the present 
work concerns with the GMA modeling of metabolic 
networks, the estimation strategies discussed in this 
work have general applicability to any kinetic models 
that can be written as X(^) = Sv(£/ C ). The creation of 
computationally efficient parameter estimation methods, 
such as the one presented here, represents an important 
step toward genome-scale kinetic modeling of cellular 
metabolism. 
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